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^-j- ' We present a general method for solving the non-linear differential equation of mono- 

tonically increasing steady-state radiation driven winds. We graphically identify all the 
singular points before transforming the momentum equation to a system of differential 
equations with all the gradients explicitly give. This permits a topological classifica- 
£->«. , tion of all singular points and to calculate the maximum and minimum mass-loss of the 

wind. We use our method to analyse for the first time the topology of the non-rotating 
frozen in ionisation m-CAK wind, with the inclusion of the finite disk correction factor 
and find up to 4 singular points, three of the x-type and one attractor-type. The only 
singular point (and solution passing through) that satisfies the boundary condition at the 
Q ' stellar surface is the standard m-CAK singular point. 

CO 

. ©2006 WTLEY-VCH Verlag GmbH & Co. KGaA, Weinheim 



^ ! 1 Introduction 

_ _ i 

Since the launch of the first satellite with a telescope on board, stellar winds from hot stars 
have been identified as a general astrophysical phenomenon. These winds are driven by the 
transfer of momentum of the radiation field to the gas by scattering of radiation in spectral 
lines. The theory of radiation driven stellar wind is the standard tool to describe the observed 
properties of the winds from these stars. Based on the Sobolev approximation, Castor, Ab- 
bott and Klein (1975, hereafter "CAK") developed an analytical hydrodynamic steady-state 
model which has been improved (or modified) later by Friend and Abbott (1986,"FA") and 
Pauldrach et al. (1986,"PPK") and general agreement with the observations could be ob- 
tained (for an extended review see Kudritzki and Puis, 2000, "KP"). The success to repro- 
duce quite well the observed winds led to the development of a new method to determine 



* Corresponding author: e-mail: michel.cure@uv.cl 



®WILEV 

I interScience* 

' ,..„■... © 2006 WILEY- VCH Verlag GmbH & Co. KGaA, Weinheim 



790 



Cure & Rial: A New Numerical Method 



extra-galactic distances, i.e. the wind momentum luminosity relation (WML,see e.g., Ku- 
dritzki, 1998, Kudritzki and Przybilla, 2003). However, in the calculations involved in the 
WML relation, the solution of the improved CAK wind (hereafter m-CAK) is not used be- 
cause the resulting velocity fields are in disagreement with the observations (KP). To avoid 
this problem usually an ad-hoc /3— field velocity profile is assumed (see KP). The main rea- 
son of this failure of the m-CAK model is very probably caused by the complex structure of 
the non-linear transcendental equation for the velocity gradient and its solution scheme. 

Due to the non-linearity in the momentum differential equation exist many solution 
branches in the integration domain. A physically reasonable solution describing the observed 
winds must start at the stellar photosphere, satisfy certain boundary condition, and finally 
reach infinity. As there is no solution branch that covers the whole integration domain, a so- 
lution must pass through a singular point that matches two different solution branches. The 
integration of the hydrodynamic differential equation describing stellar winds is therefore 
inevitably related to the existence of singular points and special numerical codes have been 
developed to identify them. In radiation driven winds two numerical algorithms to solve 
the wind's momentum equation are used. The first one uses a first guess of the location of 
the singular point (see e.g., FA or PPK) and then applies a root-finding routine to calculate 
the value of the velocity gradient. Once the velocity gradient is known, a standard routine 
(e.g. Runge-Kutta) is employee to integrate up and downstream. When the velocity field is 
attained, the lower boundary condition (see below) is calculated and the whole process is 
iterated until convergence. The main disadvantage of this method is that a initial location 
of the singular point must be known in order to integrate the equation. The second method 
(Nobili & Turolla, 1988) is a generalisation of the Henyey method (Henyey et al. 1964) to 
handle singular points. Here the nonlinear differential equation is linearised using a finite 
difference scheme. The singular point regularity condition is treated as a movable boundary 
condition (see Nobili & Turolla for details). This algorithm needs a trial velocity profile as 
initial solution and uses a Newton-Raphson method for convergence (see Krticka 2003). The 
advantage of the latter in comparison with the first one, is that no special routine is needed 
to find either the singular point location or the velocity gradient. However, the main disad- 
vantage is that the trial function must be "close" to the solution, this is important especially 
if more than one solution exists (see Cure 2004). 

The investigation of the mathematical features of nonlinear hydrodynamical equations of 
stellar winds is not new. After Parker (1960) developed the theory of the solar wind, Carovil- 
lano & King (1964) studied all the possible mathematical solutions of Parker's equation and 
clearly showed that no physical relevant solution where neglected in Parker's original analy- 
sis, i.e., only one stationary-state solution is physically acceptable. In radiation driven stellar 
winds, Bjorkman (1995) performed for the first time a detailed topological study of the orig- 
inal CAK non-linear hydrodynamic equation. He found a total of five singular points and 
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showed that the only outflow solution that satisfy the boundary condition of zero gas pres- 
sure in infinity is the CAK solution. Another step in the understanding of all possible solu- 
tions of the original CAK differential equation has been carried on by Cure & Rial (2004). 
They analysed the topology of the CAK equation with the inclusion of the centrifugal force 
term due to the star's rotational speed and also took into account the effects of changes in 
the ionisation of the wind (line-force parameter S). Indeed, another singular point could be 
found this way, but the corresponding topology is of the attractor-type and hence physically 
meaningless for outflowing winds. The work by Cure & Rial (2004) therefore confirmed the 
CAK solution being the only physical solution satisfying the lower boundary condition (see 
below) and reaching infinity, if the stellar rotation is taken into account. 

As just mentioned, the topology of the CAK equations has been analysed in quite some 
detail and several authors independently confirmed the CAK solution being the only physical 
solution. However, for the improved m-CAK model, i.e. CAK plus the finite disk correction 
term, such a detailed analysis has not yet been performed. In particular in view of the recent 
result of Cure (2004), who indeed proved the existence of additional singular points (with the 
corresponding solutions passing through) in high rotational radiation driven m-CAK winds, 
analysing in detail the topology of the m-CAK solutions appears to be highly required. In 
addition, it is crucial to understand the solution topology of the standard m-CAK model if 
one wants to incorporate other physical processes into the theory. 

In this article we present a new general numerical method that allows to solve the wind 
momentum equation in a simple and straightforward way by transforming the equations into 
a system of ordinary differential equations. A numerical (monotonicaly increasing) wind 
solution for this system can be attained using standard numerical algorithms. Then, we de- 
rive a simple condition that allows to classify the topology structure of the singular points 
in the integration domain and finally we apply this method to perform for the first time a 
topological analysis of the non-rotating frozen in ionisation m-CAK model. 

The structure of the article is the following. In section|2]we briefly review the non-linear 
differential equation for the momentum in radiation driven theory. In section [3] after intro- 
ducing a coordinate transformation, we present a graphical method to find the number and 
location(s) of the singular point(s). We introduce in section H a new physical meaningless 
independent variable to transform the wind momentum equation to a system of differential 
equations where all the derivatives (with respect to this new variable) can be obtained analyt- 
ically. In section|5]we linearise the equations in the neighbourhood of the singular point(s). 
The eigenvalues of the matrix of the linearised system give the classification topology of the 
singular point. Section[6]describes the iteration scheme to integrate the system of differential 
equations together with a lower boundary condition. In section (|7]l we apply our method by 
analysing the non-rotating frozen in ionisation radiation driven wind, including the finite 
disk correction factor (m-CAK model). 
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2 The non-linear differential equation 

The standard stationary model for line-driven stellar winds treats an one-component iso- 
thermal radial flow, ignoring the influence of heat conduction, viscosity and magnetic fields 
(see e.g., FA). For a star with mass M, radius i?*, effective temperature T e f / and luminosity 
L, the momentum equation including the centrifugal force term due to star's rotation reads: 

dv I dp GM(l-T) , u$(r) M , 

v- r = — 4 + — + g Hne (p,v',n E ) (1) 

dr p dr r z r 

where v is the fluid velocity, v' — dv/ dr is the velocity gradient, p is the mass density, p is 
the fluid pressure, v^ — v rot R*/r, where v rot is the star's rotational speed at the equator, 
r represents the ratio of the radiative acceleration due to the continuous flux mean opacity, 
er e , relative to the gravitational acceleration, T = cr e L/AircGM and the last term, g hne 7 
represents the acceleration due to the lines. 



The standard parameterisation of the line-force (Abbott, 1982) is given by 

6 

(2) 



9 Une = y 2 fn(r,v,v') (r*v v') a ^ 



tie 



W(r) 

W(r) is the dilution factor and /r>(r, v, v') is the finite disk correction factor. The line force 
parameters are: a, 5 and k (the latter is incorporated in the constant C), typical values 
of these parameter from LTE and non-LTE calculations are summarised in Lamers and 
Cassinelli (1999,"LC", chapter 8). The constant C represents the eigenvalue of the problem 
(see below) and is related to the mass loss rate (M) by: 
( 4tt \ a 

C = TGMk r- , (3) 

\(7 E VthM J 

here v t h — (2kboiT '/m^r) 1 / 2 is the hydrogen thermal speed, tie is the electron number 
density in units of 10~ n cTO~ 3 (Abbott 1982). The meaning of all other quantities is the 
standard one (see e.g., LC). 

Througout this paper, we use the original m-CAK description of the line-force. In his 
topological study of the CAK model, Bjorkman (1995) used instead the absolute value of 
the velocity gradient, to allow for the possibility of non-monotonic velocity laws. In Cure & 
Rial (2004) and in this work we focus on physically meaningful wind solutions, i.e. steady- 
state monotonically increasing outflows. 

The corresponding continuity equation reads: 

47rr 2 p v = M (4) 

Thus, replacing the density p from eq. (01, the m-CAK momentum eq. (Q~|) with the line 
force given by eq. (0 can be expressed formally as: 

F(r,v,v')=0 (5) 
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3 Singular Point Location 

In this section we describe a method to obtain the locations of the singular points. After a 
transformation of coordinates that takes advantage of the functional properties of the cor- 
rection factor, two functions are obtained. Plotting this functions in a 2-dimensional phase- 
space, after assuming a value for the constant C, the intersection points of these two curves 
give the exact locations of the singular points. 



3.1 Coordinate Transformation 

We define u, w, and w' as follows: 

u = -R*/r, (6) 

w = v/a, (7) 

w = dw/du, (8) 

where a is the isothermal sound speed. Using these new coordinates, the momentum equa- 
tion (0 transforms to: 

F(u, w, w') = (9) 
The full equation including © for the rotating m-CAK wind is given in Cure (2004, eq. [7]). 

3.2 Singular Point Conditions 

As mentioned previously, there is no solution branch that covers the entire integration do- 
main (i?* < r < oo). A physical solution must cross through a singular point, that connects 
the two solution branches. The singular point condition reads: 

JL-F(u,w,w') =F w ,(u iW ,w') = (10) 
aw' 

Hereafter we use subindex to denote partial derivatives. In addition to the singularity condi- 
tion the regularity condition has to be fulfilled, which is defined by: 

—F(u,w,w') = FJu,w,w') +w'F w (u,w,w') = (11) 
du 

3.3 Logarithmic Coordinates 

As the velocity field grows nearly exponentially near the photosphere, we transform w and 
w' to logarithmic variables. Compared to standard methods this gives a better resolution 
close to the photosphere which is important for the line formation (Venero et al. 2003). 
We define logarithmic variables: 

rj = ln(2w'/w), (12) 
C = ln(W), (13) 
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and obtain for the momentum equation (0 

F(u, V ,()=0. (14) 
The eigenvalue C is re-scaled in these coordinates, i.e., 

C = C (a 2 R*)( a -V (15) 

The singularity condition w'F w i (it, w, w') and the regularity condition F u (u, w, w')+w' F w (u, w, 
transform to: 

F v (u,r),0+F ( (u,r),0 = (16) 

F u (u, v,0 + y ( F ^ u > *>> - F ^ 0) = (17) 

Solving equations ( TBI . ([TBT l and (TTTb simultaneously we obtain the location of the criti- 
cal point u c and T) c , £ c . In radiation driven winds this is possible only if we know the value 
of the eigenvalue C, which is fixed by a boundary condition at the stellar surface. However, 
assuming C as a variable, we can find the range of values for C that allows the singular 
points to exist. Thus, we have four unknowns, u c , r) c , ( c and C and only three equations 
(1141 16l andlT7l> which we solve using the implicit function theorem. 



3.4 A Graphical Method for Locating the Singular Point(s) 

Fortunately, thanks to the line-acceleration term, it is possible to find solutions for £ in terms 
of u, T) and C, directly from the singularity condition. This feature represents the keystone 
of the methodology we develop to obtain the location of the singular points. This is also 
extendable for any additional term in the radiation driven differential equation, e.g., magnetic 
fields (Friend & McGregor, 1984) or magnetically channeled winds (Cure & Cidale, 2002; 
Owocki & ud-Doula 2004) as long as these terms do not explicitly depend on the velocity 
gradient w'. 

After solving £ — ((u, i], C) from the singularity condition and replacing it in equations 
(O and (flTl l. we obtain the following two functions: 

R(u,r},C) =0 (18) 
H(u,r),C)=0 (19) 

where (fTSl for the rotating m-CAK wind is given in Cure (2004, eq. [18]). 

Once a value of the eigenvalue C is assumed, the intersection of the above defined func- 
tions in the plane u, r\ gives the location of the singular points (see below). With the value of 
u c , we then can calculate r\ c and ( c and proceed to integrate up and downstream. 
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4 The Equivalent System of Differential Equations 

In this section, we introduce a new physical meaningless independent variable and show 
that the momentum equation transforms to a system of differential equations with all the 
derivatives explicitly given. 

Combining the definitions of r\ and C (equations [T2l and [T3l respectively) we obtain the 
following relation between its derivatives: 

f=e* + d /. (20) 
du du 

Differentiating F(u, r], C), we get: 

dF = F u du + F. n dr] + F c d( = (21) 
and using (l20l > in (fJTJ we obtain: 

dF = (F u - e"F„) du + (F„ + F c ) d( = (22) 

We now introduce a new physically meaningless independent variable, t, defined implicitly 
by: 

du= (F n + F c )dt. (23) 

With this new independent variable, equation ( TBI ) is equivalent to the following system 
of ordinary differential equations: 

^=X(u,r,,0, (24) 

^=Y(u, V ,0, (25) 

§=ZM,C)- (26) 
where X, Y and Z are defined as: 

X = F V + F (27) 

Y = -F u -e r >F ( , (28) 

Z=-F u + e^F v . (29) 

Any solution of this system of differential equations is also a solution of the original 
momentum equation (TPfl) . since if any initial condition (uq, i]o, Co) satisfies F(uo, Co) = 
0, then any solution of H24H261) verifies that F(u(t), T](t), ((t)) = 0. 

5 Topology of the Singular Points 

In this section we describe the steps involved in the topological classification of any singular 
point. 
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5.1 Linearisation 

All the critical points of the system (1241126b satisfy simultaneously F = and X = 0, 
Y = 0, Z = 0. Because the last three equations are non independent between them, it is 
sufficient to consider only two of them. 

Given (u c , r] c , £ c ), i.e. a singular point, we study the behaviour of the solutions near that 
point by considering small perturbations in its neighbourhood, thus: 

Su = X u 5u + X v 5t] + X c 8C, + o(Su,Srj,S() (30) 

Sr) = Y u Su + Y r ,Sr] + Y c 5( + o(Su,8r],8() (31) 

S( = Z u 5u + Z, 1 5r, + Z Z 5C, + o(5u,Sr],8() (32) 

^From equation ( TT4"1 >, we know Su, Sr], SC, satisfy the following constraint: 

F u Su + F n Sr] + 5( = o(5u,Sr),5() , (33) 

and therefore we can reduce this system d30l - l32l to a 2-dimensional system, i.e.: 

Su = (X u - X V FJF V ) Su + (X c - X V F C /F ri ) 6( + o (Su, Sr], 5Q (34) 

5( = (Z u - Z V F U /F V ) 6u + (Z c - Z V FJF V ) 5( + o (Su, Sr], SQ (35) 

5.2 The Matrix B, Eigenvalues and Eigenvectors 

Maintaining this expansion to first order in Su, 5( and writing the last two equations in 
matrix form, we recognise the matrix B, defined as: 

X u — X n F u / F n Xq — X^F^/ Frj 



B 



Z u — Z ri F u /F ri Zq — Z^F^/F,) 



(36) 



Because (u c , r\ c , ( c ) are defined at a singular point, we can use X = 0, Y = and Z = in 
(l36*l l. thereby obtaining: 

X u — e^Xr, Xq + X v 
Z u -e*Z n Z^ + Zr, 



B 



(37) 



(« C ,?J C ,C C ) 

The matrix B contains the information of the structure of the topology of the singular 
points. If det (B) < 0, the eigenvalues of B are real and have opposite signs. In this case, 
(u c , 7] c , (c) define a saddle type singular point (Hartman-Groessman theorem, Aman 1990), 
while the eigenvectors determine the tangent vectors in the directions of stable and unstable 
manifolds. 



6 Method of Integration 

If we start to integrate from the singular point, the values of the gradients, u, r) and ( are 
zero, because X = 0, Y = and Z = at the singular point, and therefore, we can not 
leave this point. Thus, we need to move in the neighbourhood of this singular point in the 
direction of the eigenvalues and from this new position start to integrate. 
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6.1 The stable manifold 

A physically stationary solution of the wind model must have a monotonically increasing 
velocity profile, this is topologically characterised by the condition £ > 0, thus we need 
to move in the direction of the eigenvector of B that satisfies this criterion. Defining the 
components of this eigenvector of B as (Su, 8(), we find that 6r) fulfils: 
5r)=-{F u 6u + F ( S() / F v . (38) 
Once we know the eigenvector in phase-space (u.i], (), we integrate up and downstream 
starting from the points q ± defined by: 

g ± = (u c T ?c ,Cc) ±e(5u,6r),60 (39) 
where < e « 1. The solution that starts at the initial point q~, at t = 0, achieves the 
photosphere, u = — 1 at t = T~ and the solution that starts at the initial point q + , also at 
t = 0, reaches infinity, u = at t = T + . 

6.2 Lower Boundary Condition 

The location of the singular point is determined by the lower boundary condition. In radia- 
tion driven winds from hot stars, the most frequently used boundary condition is a constraint 

on the optical depth integral: 

f 00 2 

/ g e p(r) dr=- (40) 
Jr., 6 

An equivalent lower boundary condition is to a given value of the density at the stellar 
surface, i.e. 

p{R*)=P* (41) 
Usually, in case the optical depth integral is used, it can be calculated only after the 
numerical integration has been carried on. We now show a calculation that permits to in- 
corporate the boundary condition as a fourth differential equation of the system given by 
equations ( |24"1|2"6| |. 

The boundary condition at the stellar surface is given by 
7o = J ^(x,j],C)dx (42) 
with 7q being a fixed value. Using equation ( f24b . we can write 7 as 70 = 7^ — 7^ where 

7o + = / <f>(u,Ti,C)X(u,ri,Qdt, (43) 
Jo 

and 

7o - = f <f>{u,r,,C)X{u,n,C)dt (44) 
Jo 

Then, we have 7^ = 7 (T ± ) where 7 (f) satisfies the following differential equation: 
7 = 0(x,77,C)X(x,77,C). (45) 

Thus, adding this last equation to our system of differential equations ( 12411261 ). we calcu- 
late simultaneosly the value of the lower boundary condition while the numerical integration 
is performed. 
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Table 1 B2 V Stellar Parameters 



R/R Q 


M/Mq 


L/Lq 


Teff/K 


4.5 


9.0 


3553. 


21000. 



Table 2 Line-force Parameters 



k a 8 



0.212 0.56 0.0 



6.3 The Iteration Cycle 

The iteration cycle for finding out the numerical solution of (fl4] i and satisfying the boundary 
condition ( |42T > is given by: 

- Guess a value of the Eigenvalue C 

- Calculate the values of (u c , Tj c , ( c , 7) by solving the system ( I24H26I I and d45b . and start 
to integrate numerically from (g , 0). 

- After the integration, compare the value of 70 with 7 (T + ) — 7 (T~) and modify C. 

- Repeat the cycle until convergence is reached. 

7 Topology of non-rotating frozen in ionisation m-CAK wind 

In this section we study the topology of the frozen in ionisation non-rotating m-CAK wind, 
using the methodology developed in previous sections. We solve equation (Q]i with the line 
force term given by eq. (O with 5 = 0. 

In order to compare our results with previous topological analysis (Bjorkman 1995, Cure 
& Rial 2004) we adopt the same B2 V star's parameters and line-force parameters used in 
these studies. TableQ]and[2]summarise these parameters. 

7.1 Graphical method for the location of critical points 

The first step to understand the topology of this non-linear differential equation is to know 
the location of the singular point(s). This is done using the functions i? and H, which are 
given in u, 77, £ coordinates by: 



R(u, V )= A' 1 [e" Jr) (-1 + a) f D (X) + (2A- 



f D w 

(46) 
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and 

H(u, V ,C)= -C + A-ie" (A+%) [§ e" + A~\A + f ) (a e" / D (A) - 2 u#,(A))] " 4 

(47) 

where A is an auxiliary function denned as: 



A = (l-a) e r ' f D {\)+2uf D {\) 
and 

A = M(u + 2e- ,J ) 

The analytic expression for the finite-disk correction factor is given by: 



MA) 



1 1 



l-(l-A) 



(l+a) 



(1 - a) A 

and f' D (X) denotes the total derivative of /d(A) with respect to A. 



(48) 



(49) 



(50) 




-0.8 -0.6 -0.4 



Fig. 1 The iso-contours R(u,r)) = (dashed-line) and iso-contours H(u,r),C) = 
(solid-lines) for different labelled values of C as function of the variables u and 77. See text 
for details. 



Fig. [T] shows the curves defined by R(u,T]) — (dashed lines) and H(u,T],C) = 
(solid lines), for different values of C. The intersection of this two curves gives the location 
of the singular points and the number of singular points depends on the assumed value of 
the eigenvalue. We see clearly that there are two families of curves for R(u, rj) = 0, which 
are independent of the assumed value of C. 

The standard critical point of the m-CAK model is located close to the stellar surface, 
u < — 1. On the other hand, the velocity is small near the photosphere, therefore 77 > 
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Table 3 Singular Points Coordinates 



Label C r c u c r\ c £ c 



A 
B 
C 
D 
E 



70 1.018 

46 1.116 

46 4.505 

46 29.41 

46 40.32 



-0.982 
-0.896 
-0.222 
-0.034 
-0.020 



+4.905 
-0.587 
+0.753 
-0.176 
+4.605 



+7.468 
+5.638 
+7.253 
+7.344 
+7.439 



and the family curve of the standard m-CAK model is the upper R(u, rj) = curve. The 
value of C for the singular point located in the upper R — curve (u < — 1, r\ > 0) is in 
agreement with a typical value of the mass loss rate of a B2 V star, M ~ 10~ 9 Mq yr~ Y 
which corresponds approximately to C ~ 70. The values of the coordinates u, r\ and £ at 
this critical point is given in table [3] labelled by A. 

The more complex case of C = 46 {M ~ 2 • 10~ 9 Mq yr _1 ) shows that the curves R = 
and H = intersects in four different locations, one in the upper R = curve (labelled by 
E) and three in the lower R = curve (labelled from B to D). Table [3] also summarises the 
coordinates of these critical points. 

An additional benefit of this graphical method is that this curves intersect only for a range 
of values of the eigenvalue C, therefore this method gives a lower (M m in) and an upper 
(M max ) limit of the mass loss rate of the wind. In this case (see FigureQ} it is approximately 
42 < C < 72, corresponding to 8.8 • 1(T 10 Mq yr~ x < M < 2.3 • 1(T 9 Mq yr~ x . 

7.2 The system of differential equations 

Once the location of a critical point is known we can integrate up and downstream from this 
point. As we stated, the non-linear differential equation (O can be transformed to the set of 
equations (|2~71-|29]|. In our case the functions X, Y and Z are: 



X(u,r ? ,C,C)= -Ce a ^{a.f D {\)+ 2 u e~* f' D (X)) + e< - \ e" (51) 

Y(u, V ,CC)= Ce a t (ae' /u(A) + 2 (e~ v + u) f' D (X)) - e (+v + Jr (52) 

Z(u,ri,C,C) = 2Ce a < (e-» + 2ti) f' D (\)-\e^ + ^ (53) 

Before we start to integrate, we determine the classification topology of each one of the 
singular points of Table [3] 
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Table 4 Singular points Topology classification 





det{B) 


Topology 


A 




X-type 


B 




X-type 


C 


+ 


Spiral 


D 




X-type 


E 




X-type 



7.3 Topology classification 

The Hartman-Groessman theorem (Aman 1990), states that a singular point topology can 
be obtained from the multiplication of all eigenvalues of the matrix B or equivalently its 
determinant. This matrix B is defined by eq. ( f37T >. In our case of study, B is: 

[f D (\)B (u,r))+f' D WBi(u,r)) + f£(\)B 2 (u,r))] (54) 

Where the matrix Bq is defined as: 

(A + 2/u) a 




,2 n 



(55) 



u (e 2,1 u + 2 (A+ I) n) Au (2 + Au) (-1 + a) 

4-f + 2e 3 ' ^ ^i 2 + 6e , ' (2 + Au) -2u (e 2T >u- (A+f) II) 
where the auxiliary function IT is defined as: 

n = (l-a + e* 1 (u-2ua)) (57) 
and the matrix B2 is: 

4 (2 + Au) I u(l + 2e*>u) ~u 2 \ 
e " I (l + 2e r 'u) 2 -u (l + 2e^u) J 

Table [4] summarises the value of the determinant of the B-matrix for all the singular 
points of Table [3] As expected the standard m-CAK singular point (A) is an X-type one (or 
saddle-type: eigenvalues are real, one negative and one positive). The other singular point 
(E) of the upper R = family is also an X-type singular point. The lower R = family 
of singular points, has two X-type singular points (B, D) and one spiral (C, spiral-unstable: 
complex conjugate and positive real part eigenvalues). 

Once the topology of each singular point is known, we integrate the system of differential 
equations in order to search for physical solutions of radiation driven winds. 

7.4 Numerical Integration 

Now we show the results of the numerical integration of the system of differential equations 
given by ([24]— [26j> together with d5Tl]53l ). 
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7.4.1 The Upper R = Family 

Fig. [2]shows the velocity profile as function of \og(r / R ast — 1). The first iteration starting 
from the singular points A and D are shown in dashed-lines. After performing the iteration 
described in section ( 16.3b . with the lower boundary condition given by eq. (00), both so- 
lutions converges after few iterations to the same solution, shown in Figure [2] (solid line). 
This solution corresponds to the standard m-CAK solution (FA or PPK), and the values 
of the location of the singular point, eigenvalue, mass loss rate and terminal velocity are 
respectively: r c = 1.021 i?* (u c = -0.979), C = 69.758, M = 9.39 - 10~ 10 Mq yr^ 1 , 
Voc = 2441 s _1 . 

2500 
2000 
1500 
1000 
500 

io- 2 lO" 1 10° 10 1 10 2 

r/R. -1 

Fig. 2 Velocity (in km s^ 1 ) versus (r/i?* — 1). The unique curve that starts at the stellar 
surface and reaches infinity is the m-CAK original solution (solid line). See text for details. 

7.4.2 The Lower R = Family 

Fig. [3] show the Velocity profile of the first iteration starting from the singular points B 
(dashed-line) and D (solid-line). Both solutions reach the neighbourhood of stellar surface 
significantly exceeding reasonable values of the velocity, i.e., V(R*) ~ 520 km s^ 1 and 
V(R*) ~ 960 km s" 1 when the iteration started at point B and D respectively. 

The iteration algorithm described in section ROl does not converge, because the location 
of the singular point is shifted beyond the integration domain, u > 0, for both singular 
starting points B and D. We therefore conclude that the solutions starting from singular 
point B or D do not satisfy the lower boundary condition given by eq. ( f40b . 

8 Conclusions 

We have developed a new method for solving monotonically increasing steady-state radia- 
tion driven winds which can be summarised as follows: i) find the location of critical points, 
ii) classify its topology, iii) transform the non-linear differential equations to a system of 
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Fig. 3 Velocity (in km s _1 ) versus (r/i?* — 1). The solution that passes through singular 
point B (dashed-line) reaches the photosphere with a too high value of the velocity. While 
the solution passing through the singular point D (solid-line) does not reach the photosphere. 
See text for details. 

ordinary differential equations with all the derivatives explicitly given, and finally iv) iterate 
numerically the system to obtain a wind solution. 

We applied this method to study for the first time the topology of the non-rotating, frozen 
in ionisation m-CAK wind, i.e. the CAK wind with the finite disk correction factor. In his 
CAK topological analysis, Bjorkman's (1995) found for monotonically increasing outflow 
solutions only one solution branch with one singular point and one physical wind solution, 
i.e. the original CAK solution. In Cure & Rial (2004) we identified two solution branches 
with two singular points and one physical wind solution. 

In this work we find two solution branches and up to four singular points, depending of the 
value of the Eigenvalue. The standard m-CAK singular point, however, is the only one that 
satisfies the boundary condition at the stellar surface. Furthermore, in the range of possible 
Eigenvalues, the physical solution corresponds to a value close to the maximum and there- 
fore to a minimum mass-loss rate. 

We conclude therefore, that for the boundary conditions considered here, the wind adopts 
the minimum possible mass-loss rate. 

Our new method can be easily applied to more general cases, e.g., rotating radiation 
driven winds or winds with magnetic fields. These topics are beyond the scope of this article 
and we leave them for future studies. 
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